################### PROJECT- JCR REPLICATION - figures_and_tables_appendix ##############################################
# This R file plots all tables and figures in Appendices A-T.
# Last updated: April 2025, by Wangyin Zhao
##############################################################################################################

##############################################################################################################
# Outline of Script
# A. Preparation
#    - Set working directory
#    - Load required packages
# B. Import Data
# C. Online Appendix A. Main models for the moderating impact on conflict outcomes
# D. Online Appendix B. Non-contiguous movement
# E. Online Appendix C. Spillover effects of natural disasters
# F. Online Appendix D. Parallel Trends Assumption - event study model
# G. Online Appendix E. Models that include nighttime lights as a control variable
# H. Online Appendix F. Aggregation at the municipality level
# I. Online Appendix G: Sensitivity to variations in barangay geographical size
# J. Online Appendix H: Restricting the sample to regions with NPA presence
# K. Online Appendix I: Using Conley standard error
# L. Online Appendix J: Randomization inference
# M. Online Appendix K: Matching rate between my measure and EM-DAT
# N. Online Appendix L: Accounting for rainfall concentration
# O. Online Appendix M: Zero-stage models based on EM-DAT
# P. Online Appendix N: Using alternative thresholds
# Q. Online Appendix O: A DID model accounting for the treatment phase
# R. Online Appendix R: Co-distribution of two conflict outcomes
# S. Online Appendix S: Count models for battle-related violence
# T. Online Appendix T: Exploring alternative temporal aggregations of battle-related violence
##############################################################################################################



#### A. Preparation ####
# Set working directory to Replication_Files folder
setwd("./Replication_Files")

# Import relevant packages
packages <- c("tidyverse", "sf", "broom", "xtable", "purrr", "fixest", "gridExtra", "fixest")
lapply(packages, library, character.only = TRUE)
########

#### B. Import data ####
# To generate those figures and table, the first step is to run the data processing files and to generate the analysis data 
# Upload dataset1_oldschool.csv
dataset1 <- read_csv("./data/analysis/dataset1_oldschool.csv")

# Upload dataset1_staggeredadoption.csv
dataset1_staggeredadoption <- read_csv("./data/analysis/dataset1_staggeredadoption.csv")

# Upload dataset2_oldschool.csv 
dataset2 <- read_csv("./data/analysis/dataset2_oldschool.csv")

# Upload dataset2_staggeredadoption.csv
dataset2_staggeredadoption <- read_csv("./data/analysis/dataset2_staggeredadoption.csv")

# Upload dataset1_monthlevel.csv
dataset1_month <- read_csv("./data/analysis/dataset1_month.csv")

# Upload dataset1_monthlevel.csv
dataset1_month_staggeredadoption <- read_csv("./data/analysis/dataset1_month_staggeredadoption.csv")

# Upload dataset_mun group
dataset_mun_0.2_1 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.2_1.csv")
dataset_mun_0.2_2 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.2_2.csv")
dataset_mun_0.3_1 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.3_1.csv")
dataset_mun_0.3_2 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.3_2.csv")
dataset_mun_0.4_1 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.4_1.csv")
dataset_mun_0.4_2 <- read_csv("./data/analysis/dataset_mun/dataset_mun_0.4_2.csv")

# The spatial data on the barangays from the Philippine Standard Geographic Code (PSGC) 
barangay <- st_read("./data/raw/Barangays/Barangays.shp")
########

#### Online Appendix A. Main models for the moderating impact on conflict outcomes ####
# Table A1
amod1 <- feols(Tchange_inf2 ~ ND_1.6_bar + Mixed_inf2 + Enclave_inf2 +
                ND_1.6_bar*Mixed_inf2 + ND_1.6_bar*Enclave_inf2 
              | ADM4_PCODE + year, vcov = ~ADM4_PCODE + year, 
              data = dataset1) 

amod2 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2
              | ADM4_PCODE + year, vcov = ~ADM4_PCODE + year, 
              data = dataset1_staggeredadoption) 

etable(amod1, amod2, tex = TRUE)

# Table A2
amod3 <- feols(npa_after1~ ND_1.6_bar + Mixed_inf2 + Enclave_inf2 +
                   ND_1.6_bar*Mixed_inf2 + ND_1.6_bar*Enclave_inf2
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2)

amod4 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_staggeredadoption)

etable(amod3, amod4, tex = TRUE)
########

#### Online Appendix B. Non-contiguous movement ####
# Table A3
bmod1 <- feols(Tchange_inf2 ~ treated*Influence2
             | ADM4_PCODE + year, vcov = "twoway",
             data = dataset1_staggeredadoption) 

bmod2 <- feols(Tchange_inf2 ~ treated*Mixed_inf2*Influence2 + 
               treated*Enclave_inf2*Influence2
             | ADM4_PCODE + year, vcov = "twoway", data = dataset1_staggeredadoption)

etable(bmod1, bmod2, tex = TRUE)

# Table A4
bmod3 <- feols(npa_after1 ~ treated*Infl2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption) 

bmod4 <- feols(npa_after1 ~ treated*Mixed_inf2*Infl2 + 
                 treated*Enclave_inf2*Infl2
               | ADM4_PCODE + time, vcov = "twoway", data = dataset2_staggeredadoption)

etable(bmod3, bmod4, tex = TRUE)
########

#### Online Appendix C. Spillover effects of natural disasters ####
# Table A5
cmodel1 <- feols(Tchange_inf2 ~ treated +
                   NDtreated_nb_d
                 | ADM4_PCODE + year, vcov = "twoway",
                 data = dataset1_staggeredadoption) 

cmodel2 <- feols(Tchange_inf2 ~ treated*Mixed_inf2 + treated*Enclave_inf2 + NDtreated_nb_d*Mixed_inf2_nb_d + NDtreated_nb_d*Enclave_inf2_nb_d
                 | ADM4_PCODE + year, vcov = "twoway",
                 data = dataset1_staggeredadoption) 

etable(cmodel1, cmodel2, tex = TRUE)

# Table A6
cmod3 <- feols(npa_after1~
                   treated + NDtreated_nb_d
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_staggeredadoption)

cmod4 <- feols(npa_after1~
                   treated*Mixed_inf2 + treated*Enclave_inf2 + NDtreated_nb_d*Mixed_inf2_nb_d + NDtreated_nb_d*Enclave_inf2_nb_d
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_staggeredadoption)

etable(cmod3, cmod4, tex = TRUE)
########

#### Online Appendix D. Parallel Trends Assumption - event study model #### 
#### Plot the shift of territorial control ####
# Step 1: Adjust Dataset
dataset1_month_adj <- dataset1_staggeredadoption %>%
  arrange(ADM4_PCODE, year) %>%
  group_by(ADM3_PCODE, ADM4_PCODE) %>%
  mutate(treated_lag1 = dplyr::lag(treated, 1),  # Lag 1 (12 months)
    treated_lead1 = dplyr::lead(treated, 1)) # Lead 1 (12 months)

# Step 2: Event Study Model
event_study_model <- feols(Tchange_inf2 ~ treated + treated_lag1 + treated_lead1 +
                             Mixed_inf2 + Enclave_inf2 + 
                             treated * Mixed_inf2 + treated * Enclave_inf2 +
                             treated_lag1 * Mixed_inf2 +
                             treated_lead1 * Mixed_inf2 +
                             treated_lag1 * Enclave_inf2 +
                             treated_lead1 * Enclave_inf2 |
                             ADM4_PCODE + year, vcov = "twoway",
                           data = dataset1_month_adj)

# Step 3: Extract Coefficients and Confidence Intervals
coef_table <- tidy(event_study_model) %>%
  mutate(conf.low_95 = estimate - 1.96 * std.error,  # 95% CI lower bound
    conf.high_95 = estimate + 1.96 * std.error, # 95% CI upper bound
    conf.low_90 = estimate - 1.645 * std.error, # 90% CI lower bound
    conf.high_90 = estimate + 1.645 * std.error)  # 90% CI upper bound
  

# Step 4: Filter and Prepare Coefficients for Plotting
interaction_coefs <- coef_table %>%
  filter(grepl("treated(_lag|_lead)?\\d*:.*inf2", term) | grepl("treated:.*inf2", term)) %>%
  filter(!is.na(estimate)) %>%
  mutate(
    term = factor(term, levels = c(
      # Lead terms for Mixed Setting
      "treated_lead1:Mixed_inf2", 
      # Base term for Mixed Setting
      "treated:Mixed_inf2",
      # Lag terms for Mixed Setting
      "treated_lag1:Mixed_inf2",
      # Lead terms for Enclave Setting
      "treated_lead1:Enclave_inf2",
      # Base term for Enclave Setting
      "treated:Enclave_inf2",
      # Lag terms for Enclave Setting
      "treated_lag1:Enclave_inf2")),
    setting = case_when(
      grepl("Mixed_inf2", term) ~ "Mixed Setting",
      grepl("Enclave_inf2", term) ~ "Enclave Setting")) %>%
  filter(!is.na(setting))  # Exclude NA settings if any

# Step 5
interaction_coefs <- interaction_coefs %>%
  mutate(
    estimate      = estimate * 100,
    conf.low_90   = conf.low_90 * 100,
    conf.high_90  = conf.high_90 * 100,
    conf.low_95   = conf.low_95 * 100,
    conf.high_95  = conf.high_95 * 100)

# Step 6
interaction_coefs <- interaction_coefs %>%
  mutate(time = case_when(
    grepl("lag1", term) ~ 1,
    grepl("lead1", term) ~ -1,
    grepl("treated\\:", term) & !grepl("lag1|lead1", term) ~ 0,
    TRUE ~ NA_real_))

# Step.7. Subset for Mixed Setting
mixed_setting_coefs <- interaction_coefs %>%
  filter(setting == "Mixed Setting")

# Step.8. Subset for Enclave Setting
enclave_setting_coefs <- interaction_coefs %>%
  filter(setting == "Enclave Setting")

# Step.9. 
figurea1 <- ggplot() +
  # Enclave panel
  geom_pointrange(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_95, ymax = conf.high_95, color = setting),
    alpha = 0.5, linewidth = 0.5) +
  geom_pointrange(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_90, ymax = conf.high_90, color = setting),
    linewidth = 1) +
  geom_line(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, color = setting),
    linewidth = 0.5, alpha = 0.8) + 
  # Mixed panel
  geom_pointrange(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_95, ymax = conf.high_95, color = setting),
    alpha = 0.5, linewidth = 0.5) +
  geom_pointrange(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_90, ymax = conf.high_90, color = setting),
    linewidth = 1) +
  geom_line(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, color = setting),
    linewidth = 0.5, alpha = 0.8) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  facet_wrap(~setting, scales = "free_y") +
  scale_x_continuous(
    breaks = c(-1, 0, 1),
    labels = c("-1", "0", "1")) +
  labs(
    x = "Years Since Natural Disasters",
    y = " % Change in a Shift in Territorial Control",
    color = "Setting") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "top") 
########

#### Plot the battle-related violence ####
# Step 1: Adjust Dataset
dataset1_month_adj <- dataset2_staggeredadoption %>%
  arrange(ADM4_PCODE, year, month) %>%
  group_by(ADM4_PCODE) %>%
  mutate(treated_lag1 = dplyr::lag(treated, 1),
         treated_lag2 = dplyr::lag(treated, 2),
         treated_lag3 = dplyr::lag(treated, 3),
         treated_lag4 = dplyr::lag(treated, 4),
         treated_lag5 = dplyr::lag(treated, 5),
         treated_lag6 = dplyr::lag(treated, 6),
         treated_lead1 = dplyr::lead(treated, 1),
         treated_lead2 = dplyr::lead(treated, 2),
         treated_lead3 = dplyr::lead(treated, 3),
         treated_lead4 = dplyr::lead(treated, 4),
         treated_lead5 = dplyr::lead(treated, 5),
         treated_lead6 = dplyr::lead(treated, 6)) 

# Step 2: Event Study Model
event_study_model <- feols(npa ~ treated + treated_lag1 + treated_lag2 + treated_lag3 + treated_lag4 + treated_lag5  +treated_lag6+
                             treated_lead1 + treated_lead2 + treated_lead3 + treated_lead4 + treated_lead5   +treated_lead6+
                             Mixed_inf2 + Enclave_inf2 +
                             treated*Mixed_inf2 + treated*Enclave_inf2 +
                             treated_lag1*Mixed_inf2 + treated_lag2*Mixed_inf2 +
                             treated_lag3*Mixed_inf2 + treated_lag4*Mixed_inf2 +
                             treated_lag5*Mixed_inf2  +treated_lag6*Mixed_inf2  +
                             treated_lead1*Mixed_inf2 + treated_lead2*Mixed_inf2 +
                             treated_lead3*Mixed_inf2 + treated_lead4*Mixed_inf2 +
                             treated_lead5*Mixed_inf2  +treated_lead6*Mixed_inf2+
                             treated_lag1*Enclave_inf2 + treated_lag2*Enclave_inf2 +
                             treated_lag3*Enclave_inf2 + treated_lag4*Enclave_inf2 +
                             treated_lag5*Enclave_inf2  + treated_lag6*Enclave_inf2+
                             treated_lead1*Enclave_inf2 + treated_lead2*Enclave_inf2 +
                             treated_lead3*Enclave_inf2 + treated_lead4*Enclave_inf2 +
                             treated_lead5*Enclave_inf2 +treated_lead6*Enclave_inf2 |
                             ADM4_PCODE + time,
                           vcov = "twoway",
                           data = dataset1_month_adj)

# Step 3: Extract Coefficients and Confidence Intervals
coef_table <- tidy(event_study_model) %>%
  mutate(conf.low_95 = estimate - 1.96 * std.error,  # 95% CI lower bound
         conf.high_95 = estimate + 1.96 * std.error, # 95% CI upper bound
         conf.low_90 = estimate - 1.645 * std.error, # 90% CI lower bound
         conf.high_90 = estimate + 1.645 * std.error)  # 90% CI upper bound


# Step 4: Filter and Prepare Coefficients for Plotting
interaction_coefs <- coef_table %>%
  filter(grepl("treated(_lag|_lead)?\\d*:.*inf2", term) | grepl("treated:.*inf2", term)) %>%
  filter(!is.na(estimate)) %>%
  mutate(term = factor(term, levels = c("treated_lead6:Mixed_inf2",
                                    "treated_lead5:Mixed_inf2", 
                                    "treated_lead4:Mixed_inf2", "treated_lead3:Mixed_inf2", 
                                    "treated_lead2:Mixed_inf2", "treated_lead1:Mixed_inf2", 
                                    "treated:Mixed_inf2",
                                    "treated_lag1:Mixed_inf2", "treated_lag2:Mixed_inf2", 
                                    "treated_lag3:Mixed_inf2", "treated_lag4:Mixed_inf2", 
                                    "treated_lag5:Mixed_inf2",  "treated_lag6:Mixed_inf2",
                                    "treated_lead5:Enclave_inf2","treated_lead6:Enclave_inf2",
                                    "treated_lead4:Enclave_inf2", "treated_lead3:Enclave_inf2", 
                                    "treated_lead2:Enclave_inf2", "treated_lead1:Enclave_inf2", 
                                    "treated:Enclave_inf2",
                                    "treated_lag1:Enclave_inf2", "treated_lag2:Enclave_inf2", 
                                    "treated_lag3:Enclave_inf2", "treated_lag4:Enclave_inf2", 
                                    "treated_lag5:Enclave_inf2", "treated_lag6:Enclave_inf2")),
    setting = case_when(
      grepl("Mixed_inf2", term) ~ "Mixed Setting",
      grepl("Enclave_inf2", term) ~ "Enclave Setting")) %>%
  filter(!is.na(setting))  # Exclude NA settings if any

# Step 5
interaction_coefs <- interaction_coefs %>%
  mutate(
    estimate      = estimate * 100,
    conf.low_90   = conf.low_90 * 100,
    conf.high_90  = conf.high_90 * 100,
    conf.low_95   = conf.low_95 * 100,
    conf.high_95  = conf.high_95 * 100)

# Step 6
interaction_coefs <- interaction_coefs %>%
  mutate(time = case_when(
    grepl("lag1", term) ~ 1,
    grepl("lag2", term) ~ 2,
    grepl("lag3", term) ~ 3,
    grepl("lag4", term) ~ 4,
    grepl("lag5", term) ~ 5,
    grepl("lag6", term) ~ 6,
    grepl("lead1", term) ~ -1,
    grepl("lead2", term) ~ -2,
    grepl("lead3", term) ~ -3,
    grepl("lead4", term) ~ -4,
    grepl("lead5", term) ~ -5,
    grepl("lead6", term) ~ -6,
    grepl("treated\\:", term) & !grepl("lag1|lead1", term) ~ 0,
    TRUE ~ NA_real_))

# Step.7. Subset for Mixed Setting
mixed_setting_coefs <- interaction_coefs %>%
  filter(setting == "Mixed Setting")

# Step.8. Subset for Enclave Setting
enclave_setting_coefs <- interaction_coefs %>%
  filter(setting == "Enclave Setting")

# Step.9. 
figurea2 <- ggplot() +
  # Enclave panel
  geom_pointrange(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_95, ymax = conf.high_95, color = setting),
    alpha = 0.5, linewidth = 0.5) +
  geom_pointrange(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_90, ymax = conf.high_90, color = setting),
    linewidth = 1) +
  geom_line(
    data = enclave_setting_coefs,
    aes(x = time, y = estimate, color = setting),
    linewidth = 0.5, alpha = 0.8) + 
  # Mixed panel
  geom_pointrange(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_95, ymax = conf.high_95, color = setting),
    alpha = 0.5, linewidth = 0.5) +
  geom_pointrange(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, ymin = conf.low_90, ymax = conf.high_90, color = setting),
    linewidth = 1) +
  geom_line(
    data = mixed_setting_coefs,
    aes(x = time, y = estimate, color = setting),
    linewidth = 0.5, alpha = 0.8) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  facet_wrap(~setting, scales = "free_y") +
  scale_x_continuous(
    breaks = c(-6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6),
    labels = c("-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) +
  labs(
    x = "Months Since Natural Disasters",
    y = " % Change in Battle-related Violence",
    color = "Setting") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "top") 
########
########

#### Online Appendix E. Models that include nighttime lights as a control variable ####
# Table A7
emod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2 + NightLight 
              | ADM4_PCODE +year, vcov = ~ADM4_PCODE + year, 
              data = dataset1_staggeredadoption) 

emod2 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2 + NightLight 
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_staggeredadoption)

etable(emod1, emod2, tex = TRUE)
########

#### Online Appendix F. Aggregation at the municipality level ####
# Table A8
fmod1 <- feols(Tchange_in2_0.2 
                ~ nd_1.6*Mixed_inf2_0.2 + nd_1.6*Enclave_inf2_0.2
                | ADM3_PCODE + year, vcov = "twoway",
                data = dataset_mun_0.2_1)

fmod2 <- feols(Tchange_in2_0.3 
               ~ nd_1.6*Mixed_inf2_0.3 + nd_1.6*Enclave_inf2_0.3
               | ADM3_PCODE + year, vcov = "twoway",
               data = dataset_mun_0.3_1)

fmod3 <- feols(Tchange_in2_0.4 
               ~ nd_1.6*Mixed_inf2_0.4 + nd_1.6*Enclave_inf2_0.4
               | ADM3_PCODE + year, vcov = "twoway",
               data = dataset_mun_0.4_1)

etable(fmod1, fmod2, fmod3, tex = TRUE)

# Table A9
fmod4 <- feols(Tchange_in2_0.2 
               ~ nd*Mixed_inf2_0.2 + nd*Enclave_inf2_0.2
               | ADM3_PCODE + year, vcov = "twoway",
               data = dataset_mun_0.2_2)

fmod5 <- feols(Tchange_in2_0.3 
               ~ nd*Mixed_inf2_0.3 + nd*Enclave_inf2_0.3
               | ADM3_PCODE + year, vcov = "twoway",
               data = dataset_mun_0.3_2)

fmod6 <- feols(Tchange_in2_0.4 
               ~ nd*Mixed_inf2_0.4 + nd*Enclave_inf2_0.4
               | ADM3_PCODE + year, vcov = "twoway",
               data = dataset_mun_0.4_2)

etable(fmod4, fmod5, fmod6, tex = TRUE)
########

#### Online Appendix G: Sensitivity to variations in barangay geographical size ####
# Table 10
phl <- control %>% st_make_valid()
phl <- phl %>% mutate(area = st_area(phl),
                      z_area = (area-mean(area))/sd(area))

quantiles <- quantile(phl$z_area, probs = seq(0, 1, 0.05), na.rm = TRUE)

quantiles_df <- data.frame(
  `Quantile` = names(quantiles),
  `Z-Score` = as.numeric(quantiles)) #make latex table

latex_table <- xtable(quantiles_df, digits = c(0, 0, 3))

# Table A11
phl <- phl %>% st_make_valid()

phl_10 <- phl %>% mutate(z_area = as.numeric(z_area)) %>%
  filter(z_area <= -0.39695999 | z_area >= 1.00235215) %>% mutate(ADM4_PCODE_b = str_replace(ADM4_PCODE, "^PH0*", ""))

phl_5 <- phl %>% mutate(z_area = as.numeric(z_area)) %>%
  filter(z_area <= -0.39878709 | z_area >= 1.87515409) %>% mutate(ADM4_PCODE_b = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset1_filtered_2010_2011_b <- dataset1_staggeredadoption %>%
  filter(!ADM4_PCODE %in% phl_10$ADM4_PCODE)

gmod1 <- feols(Tchange_inf2 ~ 
                  treated*Mixed_inf2 + treated*Enclave_inf2
                | ADM4_PCODE + year, vcov = "twoway",
                data = dataset1_filtered_2010_2011_b) 

dataset1_filtered_2010_2011_b <- dataset1_staggeredadoption %>%
  filter(!ADM4_PCODE %in% phl_5$ADM4_PCODE)


gmod2 <- feols(Tchange_inf2 ~ 
                  treated*Mixed_inf2 + treated*Enclave_inf2
                | ADM4_PCODE + year, vcov = "twoway",
                data = dataset1_filtered_2010_2011_b) 
etable(gmod1, gmod2, tex = TRUE)


# Table A12
dataset2_filtered_b <- dataset2_staggeredadoption %>%
  filter(!ADM4_PCODE %in% phl_10$ADM4_PCODE_b)


gmod3 <- feols(npa_after1~ 
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_filtered_b)


dataset2_filtered_b <- dataset2_staggeredadoption %>%
  filter(!ADM4_PCODE %in% phl_5$ADM4_PCODE_b)

gmod4 <- feols(npa_after1~ 
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, vcov = "twoway",
                 data = dataset2_filtered_b)
etable(gmod3, gmod4, tex = TRUE)
########
 
#### Online Appendix H: Restricting the sample to regions with NPA presence ####
# Model1 in Table A13
# Step.1. Identify the adm3 with the NPA presence
dataset1_f <- dataset1 %>% group_by(ADM3_PCODE) %>%
   summarise(Influence2 = sum(Influence2, na.rm = TRUE))
 
names <- unique(dataset1_f$ADM3_PCODE[dataset1_f$Influence2 >0])
 
# Step.2. Choose the adm3 with the NPA presence only 
dataset1_t <- dataset1_staggeredadoption %>% filter(ADM3_PCODE %in% names)
 
mod1 <- feols(Tchange_inf2 ~ 
                 treated*Mixed_inf2 + treated*Enclave_inf2
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_t) 
 
# Model2 in Table A13
# Step.1. Add the ADM2
map_bar <- barangay %>% st_make_valid()
map_bar <- map_bar %>% dplyr::select(ADM2_PCODE,ADM3_PCODE) %>% st_drop_geometry() 
map_bar <- distinct(map_bar)
 
dataset_p <- dataset1_staggeredadoption %>% left_join(map_bar)
 
dataset_p_f <- dataset_p %>% group_by(ADM2_PCODE) %>%
   summarise(Influence2 = sum(Influence2, na.rm = TRUE))
names <- unique(dataset_p_f$ADM2_PCODE[dataset_p_f$Influence2 >0])
 
datasetp_t <- dataset_p %>% filter(ADM2_PCODE %in% names)
 
 
mod2 <- feols(Tchange_inf2 ~ 
                 treated*Mixed_inf2 + treated*Enclave_inf2
               | ADM4_PCODE + year, vcov = "twoway",
               data = datasetp_t) 

etable(mod1, mod2,
       tex = TRUE,           
       arraystretch = 0.6, title = "The Effect of Natural Disasters on the Shift of Territorial Control")
########

#### Online Appendix I: Using Conley standard error ####
# Table A14
dataset1_staggeredadoption <- dataset1_staggeredadoption %>% rename(lat = LATITUDE,
                                                                    lon = LONGITUDE) #adjust the name for conley

Imod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE +year, conley(10), 
                 data = dataset1_staggeredadoption) 
  
Imod2 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE +year, conley(20), 
                 data = dataset1_staggeredadoption) 
  
Imod3 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE +year, conley(50), 
                 data = dataset1_staggeredadoption) 
  
Imod4 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE +year, conley(100), 
                 data = dataset1_staggeredadoption) 
  
Imod5 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE +year, conley(200), 
                 data = dataset1_staggeredadoption) 
  
etable(Imod1, Imod2, Imod3, Imod4, Imod5, tex = TRUE)

# Table A15
dataset2_staggeredadoption <- dataset2_staggeredadoption %>% rename(lat = LATITUDE,
                                                                    lon = LONGITUDE) #adjust the name for conley
Imod6 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
               treated*Mixed_inf2 + treated*Enclave_inf2
              | ADM4_PCODE + time, conley(10),
               data = dataset2_staggeredadoption)
  
Imod7 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
               treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, conley(20),
               data = dataset2_staggeredadoption)
  
Imod8 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, conley(50),
                 data = dataset2_staggeredadoption)
  
Imod9 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, conley(100),
                 data = dataset2_staggeredadoption)
  
Imod10 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                   treated*Mixed_inf2 + treated*Enclave_inf2
                 | ADM4_PCODE + time, conley(200),
                 data = dataset2_staggeredadoption)

etable(Imod6, Imod7, Imod8, Imod9, Imod10, tex = TRUE)
########

#### Online Appendix J: Randomization inference #### 
#### Figure A3 ####
# Step.1. Define reshuffling function at the year level
reshuffle_treatment_year <- function(data, treatment_col, year_col, unit_col) {
  data %>%
    group_by(!!sym(unit_col)) %>%
    mutate(!!sym(treatment_col) := sample(!!sym(treatment_col), replace = FALSE)) %>%
    ungroup()} 

# Step.2. Estimate the original model
mod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2 |
                ADM4_PCODE + year, vcov = "twoway", data = dataset1_staggeredadoption)
original_effect1 <- coef(mod1)["treated:Mixed_inf2"] 

# Run 1000 reshuffles and estimate placebo effects
set.seed(123)  # For reproducibility
placebo_effects1 <- map_dbl(1:1000, function(i) {
  reshuffled_data <- reshuffle_treatment_year(dataset1_staggeredadoption, "treated", "year", "ADM4_PCODE")
  mod <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2  
                  |
                 ADM4_PCODE + year,  vcov = ~ADM4_PCODE + year, data = reshuffled_data)
  coef(mod)["treated:Mixed_inf2"]})

# Step.3. Compare original effect to placebo distribution
p_value <- mean(placebo_effects1 >= original_effect1)

# Step.4. Plot placebo effects distribution
hist(placebo_effects1, breaks = 30, main = "Placebo Effects Distribution",
     xlab = "Placebo Treatment Effects", col = "lightblue", border = "white")
abline(v = original_effect1, col = "red", lwd = 2, lty = 2)  # Add original effect line
########

#### Figure A4 ####
# Step.1. Define reshuffling function at the year level
reshuffle_treatment_year <- function(data, treatment_col, year_col, unit_col) {
  data %>%
    group_by(!!sym(unit_col)) %>%
    mutate(!!sym(treatment_col) := sample(!!sym(treatment_col), replace = FALSE)) %>%
    ungroup()} 

# Step.2. Estimate the original model
mod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2  
                 |
                ADM4_PCODE + year, vcov = "twoway", data = dataset1_staggeredadoption)
original_effect1 <- coef(mod1)["treated:Enclave_inf2"] 

# Run 1000 reshuffles and estimate placebo effects
set.seed(130)  # For reproducibility
placebo_effects1 <- map_dbl(1:1000, function(i) {
  reshuffled_data <- reshuffle_treatment_year(dataset1_staggeredadoption, "treated", "year", "ADM4_PCODE")
  mod <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2 
                  |
                 ADM4_PCODE + year,  vcov = ~ADM4_PCODE + year, data = reshuffled_data)
  coef(mod)["treated:Enclave_inf2"]})

# Step.3. Compare original effect to placebo distribution
p_value <- mean(placebo_effects1 >= original_effect1)

# Step.4. Plot placebo effects distribution
hist(placebo_effects1, breaks = 30, main = "Placebo Effects Distribution",
     xlab = "Placebo Treatment Effects", col = "lightblue", border = "white")
abline(v = original_effect1, col = "red", lwd = 2, lty = 2)  # Add original effect line
########

#### Figure A5 ####
# Step.1. Define reshuffling function at the month (time) level
reshuffle_treatment_time <- function(data, treatment_col, time_col, unit_col) {
  data %>%
    group_by(!!sym(unit_col)) %>%
    mutate(!!sym(treatment_col) := sample(!!sym(treatment_col), replace = FALSE)) %>%
    ungroup()}

# Step.2. Estimate the original model
mod2 <- feols(npa_after1 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2 |
                ADM4_PCODE + time, vcov = "twoway", data = dataset2_staggeredadoption)

original_effect <- coef(mod2)["treated:Mixed_inf2"]

# Step.3. Run 1000 reshuffles and estimate placebo effects
set.seed(123)  # For reproducibility
placebo_effects <- map_dbl(1:1000, function(i) {
  reshuffled_data <- reshuffle_treatment_time(dataset2_staggeredadoption, "treated", "time", "ADM4_PCODE")
  mod <- feols(npa_after1 ~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2 |
                 ADM4_PCODE + time, vcov = "twoway", data = reshuffled_data)
  coef(mod)["treated:Mixed_inf2"]})

# Step.4. Compare original effect to placebo distribution
p_value <- mean(placebo_effects >= original_effect)

# Step.5. Plot placebo effects distribution
hist(placebo_effects, breaks = 30,
     xlab = "Placebo Treatment Effects", col = "lightblue", border = "white", main = "",
     xlim = c(-0.01, 0.03))

abline(v = original_effect, col = "red", lwd = 2, lty = 2)  # Add original effect line

text(x = original_effect, 
     y = max(hist(placebo_effects, plot = FALSE)$counts) * 0.8,  # Dynamic y position, 80% of max height
     labels = paste0("Original Effect: ", round(original_effect, 3)), 
     col = "red", pos = 4)  # pos = 4 positions text to the right of the line
########

#### Figure A6 ####
# Step.1. Define reshuffling function at the month (time) level
reshuffle_treatment_time <- function(data, treatment_col, time_col, unit_col) {
  data %>%
    group_by(!!sym(unit_col)) %>%
    mutate(!!sym(treatment_col) := sample(!!sym(treatment_col), replace = FALSE)) %>%
    ungroup()}

# Step.2. Estimate the original model
mod2 <- feols(npa_after1 ~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2 |
                ADM4_PCODE + time, vcov = "twoway", data = dataset2_staggeredadoption)

original_effect <- coef(mod2)["treated:Enclave_inf2"]

# Step.3. Run 1000 reshuffles and estimate placebo effects
set.seed(125)  # For reproducibility
placebo_effects <- map_dbl(1:1000, function(i) {
  reshuffled_data <- reshuffle_treatment_time(dataset2_staggeredadoption, "treated", "time", "ADM4_PCODE")
  mod <- feols(npa_after1 ~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2 |
                 ADM4_PCODE + time, vcov = "twoway", data = reshuffled_data)
  coef(mod)["treated:Enclave_inf2"]})

# Step.4. Compare original effect to placebo distribution
p_value <- mean(placebo_effects >= original_effect)

# Step.5. Plot placebo effects distribution
hist(placebo_effects, breaks = 30,
     xlab = "Placebo Treatment Effects", col = "lightblue", border = "white", main = "")

abline(v = original_effect, col = "red", lwd = 2, lty = 2)  # Add original effect line

text(x = original_effect, 
     y = max(hist(placebo_effects, plot = FALSE)$counts) * 0.8,  # Dynamic y position, 80% of max height
     labels = paste0("Original Effect: ", round(original_effect, 3)), 
     col = "red", pos = 4)  # pos = 4 positions text to the right of the line
########
#########

##### Online Appendix K: Matching rate between my measure and EM-DAT ####
dataset_test <- st_read("data/precision_test.csv") %>%
  separate(col = code, 
           into = c("ADM4_PCODE", "year", "month"), 
           sep = "_")

dataset_test <- distinct(dataset_test)

nd <- st_read("/Users/zhaowangyin/Desktop/rainfall deviation_bar.csv") %>% 
  mutate(z_value = as.numeric(z_value),
         ND_1.6_bar = if_else(z_value >= 1.6, 1, 0)) %>% 
  rename(z_value_bar = `z_value`) %>%
  na.omit()

nd <- nd %>% filter(year > 2011 & year < 2015)

dataset_test <- dataset_test %>% left_join(nd) %>% mutate(fit_1.6 = if_else(ND == ND_1.6_bar, 1, 0)) %>%
  filter(year > 2011 & year < 2015)

#to check the comparison between em-dat and mine
## prop.table(table(dataset_test$fit_1.6))

#to summarise natural disasters across year 
prop_table <- prop.table(table(dataset_test$fit_1.6, dataset_test$year), margin = 2)

rounded_prop_table <- round(prop_table, digits = 4)

row_sums <- rowSums(table(dataset_test$fit_1.6, dataset_test$year))

col_sums <- colSums(table(dataset_test$fit_1.6, dataset_test$year))

total_count <- sum(row_sums)

row_proportions <- row_sums / total_count

row_proportions <- round(row_proportions, 4)

prop_table_with_sums <- cbind(rounded_prop_table, Row_Sum = row_proportions)

colnames(prop_table_with_sums)[ncol(prop_table_with_sums)] <- "Total"

rownames(prop_table_with_sums) <- c("Unmatched", "Matched")

kable(prop_table_with_sums, format = "latex", align = "c", 
      caption = "The Matching Rate between my Measure and EM-DAT",
      col.names = c("2012", "2013", "2014", "Total"),
      digits = 4) 
########

#### Online Appendix L: Accounting for rainfall concentration ####
# Table A17
lmod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2
               | ADM4_PCODE +year, vcov = ~ADM4_PCODE + year, 
               data = dataset1_month_staggeredadoption) 

lmod2 <- feols(Tchange_inf2 ~ treated + rci + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2*rci + treated*Enclave_inf2*rci
               | ADM4_PCODE +year, vcov = ~ADM4_PCODE + year, 
               data = dataset1_month_staggeredadoption) 
etable(lmod1, lmod2, tex = TRUE)

# Table A18
lmod3 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                treated*Mixed_inf2 + treated*Enclave_inf2
              | ADM4_PCODE + time, vcov = "twoway",
              data = dataset2_staggeredadoption)

lmod4 <- feols(npa_after1~ 
                treated*Mixed_inf2*rci + treated*Enclave_inf2*rci
              | ADM4_PCODE + time, vcov = "twoway",
              data = dataset2_staggeredadoption)
etable(lmod3, lmod4, tex = TRUE)
########

#### Online Appendix M: Zero-stage models based on EM-DAT ####
# Table A19
stage1_1 <- feols(nd_em_dat ~ z_value_bar  + rci  + I(z_value_bar > 0.7) +  I(z_value_bar > 1.6) | 
                  ADM4_PCODE + month , 
                vcov = "twoway",
                data = dataset1_month)

dataset1_month <- dataset1_month %>%
  mutate(Predicted.Disaster = predict(stage1_1, newdata = dataset1_month))

mmod1 <- feols(Tchange_inf2 ~ Predicted.Disaster + Mixed_inf2 + Enclave_inf2 + 
               Predicted.Disaster*Mixed_inf2 + Predicted.Disaster*Enclave_inf2 | ADM4_PCODE + year, 
             vcov = ~ADM4_PCODE + year, 
             data = dataset1_month)

etable(stage1_1, mmod1, tex = TRUE)

# Table A20
stage1_2 <- feglm(nd_em_dat ~ z_value_bar + rci + wet_season + 
                  I(z_value_bar > 0.7) + I(z_value_bar > 1.6) | 
                  ADM4_PCODE + month,
                vcov = "twoway",
                family = binomial(link = "logit"), 
                data = dataset1_month)

dataset1_month <- dataset1_month %>%
  mutate(Predicted.Disaster.Dummy = ifelse(predict(stage1_2, newdata = dataset1_month, type = "response") > 0.5, 1, 0))

mmod2 <- feols(Tchange_inf2 ~ Predicted.Disaster.Dummy + Mixed_inf2 + Enclave_inf2 + 
                 Predicted.Disaster.Dummy*Mixed_inf2 + Predicted.Disaster.Dummy*Enclave_inf2 | ADM4_PCODE + year, 
             vcov = ~ADM4_PCODE + year, 
             data = dataset1_month)

etable(stage1_2, mmod2, tex = TRUE)

# Table A21
stage1_3 <- feols(nd_em_dat ~ z_value_bar  + rci + I(z_value_bar > 0.7) + I(z_value_bar > 1.6) | 
                  ADM4_PCODE + month, 
                vcov = "twoway",
                data = dataset2)

dataset2 <- dataset2 %>%
  mutate(Predicted.Disaster = predict(stage1_3, newdata = dataset2))


mmodel3 <- feols(npa_after1 ~ Predicted.Disaster + Mixed_inf2 + Enclave_inf2 + 
                    Predicted.Disaster*Mixed_inf2 + Predicted.Disaster*Enclave_inf2  | ADM4_PCODE + time, 
                  vcov = "twoway", 
                  data = dataset2)

etable(stage1_3, mmodel3, tex = TRUE)

# Table A22

stage1_4 <- feglm(nd_em_dat ~ z_value_bar + rci + 
                    I(z_value_bar > 0.7) + I(z_value_bar > 1.6) | 
                    ADM4_PCODE + month,
                  vcov = "twoway",
                  family = binomial(link = "logit"), 
                  data = dataset2)

dataset2 <- dataset2 %>%
  mutate(Predicted.Disaster.Dummy = ifelse(predict(stage1_4, newdata = dataset2, type = "response") > 0.5, 1, 0))

mmod4 <- feols(npa_after1 ~ Predicted.Disaster.Dummy + Mixed_inf2 + Enclave_inf2 + 
                    Predicted.Disaster.Dummy*Mixed_inf2 + Predicted.Disaster.Dummy*Enclave_inf2  | ADM4_PCODE + time, 
                  vcov = "twoway", 
                  data = dataset2)
etable(stage1_4, mmod4, tex = TRUE)
########

#### Online Appendix N: Using alternative thresholds ####
#### Figure A7 ####
# z-value = 1.5
dataset1 <- dataset1 %>%
  mutate(year = as.numeric(year),
         numeric_time = (year - 2012))

# Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.5_bar == 1], na.rm = TRUE),  # First disaster
    treated1.5 = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

# Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.5_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.5)
# to get the data of output_f, run the data_processing1 first
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())

dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)

main1 <- feols(Tchange_inf2 ~ treated1.5*Mixed_inf2 +
                 treated1.5*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# z-value = 1.6
dataset1 <- dataset1 %>%
  mutate(year = as.numeric(year),
         numeric_time = (year - 2012))

## Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.6_bar == 1], na.rm = TRUE),  # First disaster
    treated = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

## Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.6_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.6)
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())

dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)

main2 <- feols(Tchange_inf2 ~ treated*Mixed_inf2 +
                 treated*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# z-value = 1.7
# Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.7_bar == 1], na.rm = TRUE),  # First disaster
    treated1.7 = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

## Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.7_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.7)
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())

dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)


main3 <- feols(Tchange_inf2 ~ treated1.7*Mixed_inf2 +
                 treated1.7*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# z-value = 1.8
# Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.8_bar == 1], na.rm = TRUE),  # First disaster
    treated1.8 = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

# Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.8_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)


test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.8)
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())

dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)


main4 <- feols(Tchange_inf2 ~ treated1.8*Mixed_inf2 +
                 treated1.8*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# z-value = 1.9
# Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.9_bar == 1], na.rm = TRUE),  # First disaster
    treated1.9 = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

# Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.9_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.9)
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())

dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)

main5 <- feols(Tchange_inf2 ~ treated1.9*Mixed_inf2 +
                 treated1.9*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# z-value = 2
# Code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_2_bar == 1], na.rm = TRUE),  # First disaster
    treated2 = ifelse(numeric_time >= first_disaster_year, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

# Identify second disaster time
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_2_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

dataset1_exlcude_second <- dataset1 %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 2)
test1 <- test %>% group_by(ADM4_PCODE) %>%summarise(n=n())


dataset1_filtered_2010_2011 <- dataset1_exlcude_second %>%
  filter(!ADM4_PCODE %in% test1$ADM4_PCODE)

main6 <- feols(Tchange_inf2 ~ treated2*Mixed_inf2 +
                 treated2*Enclave_inf2 
               | ADM4_PCODE + year, vcov = "twoway",
               data = dataset1_filtered_2010_2011)

# Begin to plot 
model_list <- list(main2, main1, main3, main4, main5, main6)
pd <- position_dodge(0.75)

# Extract coefficients and CIs using the broom package
coefficients_df_main <- bind_rows(lapply(model_list, tidy))
coefficients_df_main1 <- coefficients_df_main[c(4,9,14,19,24,29),]

# Create the coefficient plot with CIs using ggplot2
figureA4 <- ggplot(coefficients_df_main1, aes(y = estimate, x = term)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Coefficient estimate (mixed setting)") +
  scale_x_discrete(labels = c("This research", "1.5", "1.7", "1.8", "1.9", "2")) +
  theme_bw(base_size = 14)  +
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6) 

# Extract coefficients and CIs using the broom package
coefficients_df_main2 <- coefficients_df_main[c(5,10,15,20,25,30),]

# Create the coefficient plot with CIs using ggplot2
figureA5 <- ggplot(coefficients_df_main2, aes(y = estimate, x = term)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Coefficient estimate (enclave setting)") +
  scale_x_discrete(labels = c("This research", "1.5", "1.7", "1.8", "1.9", "2")) +
  theme_bw(base_size = 14)  +
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6) 

figure7 <- grid.arrange(
  figureA4, figureA5)
########

#### Figure A8 ####
dataset2_f <- dataset2 %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month)  # Convert to numeric representation

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_1.5_bar == 1], na.rm = TRUE),  # First disaster
    treated1.5 = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.5_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.5)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main1 <- feols(npa_after1 ~ treated1.5 + Mixed_inf2 + Enclave_inf2 +
                 treated1.5*Mixed_inf2 + treated1.5*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

dataset2_f <- dataset2_f %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month)  

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_1.6_bar == 1], na.rm = TRUE),  # First disaster
    treated = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.6_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.6)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main2 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                 treated*Mixed_inf2 + treated*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

## Code the first hit 

dataset2_f <- dataset2_f %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month)  # Convert to numeric representation


dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_1.7_bar == 1], na.rm = TRUE),  # First disaster
    treated1.7 = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.7_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()


# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.7)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main3 <- feols(npa_after1~ treated1.7 + Mixed_inf2 + Enclave_inf2 +
                 treated1.7*Mixed_inf2 + treated1.7*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

## Code the first hit 

dataset2_f <- dataset2_f %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month)  # Convert to numeric representation

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_1.8_bar == 1], na.rm = TRUE),  # First disaster
    treated1.8 = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.8_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.8)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main4 <- feols(npa_after1~ treated1.8 + Mixed_inf2 + Enclave_inf2 +
                 treated1.8*Mixed_inf2 + treated1.8*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

## Code the first hit 

dataset2_f <- dataset2 %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month  # Convert to numeric representation
  )

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_1.9_bar == 1], na.rm = TRUE),  # First disaster
    treated1.9 = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_1.9_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 1.9)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main5 <- feols(npa_after1~ treated1.9 + Mixed_inf2 + Enclave_inf2 +
                 treated1.9*Mixed_inf2 + treated1.9*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

## Code the first hit 

dataset2_f <- dataset2_f %>%
  mutate(year = as.numeric(year),
         month = as.numeric(month),
         numeric_time = (year - 2012) * 12 + month  # Convert to numeric representation
  )

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year_month = min(numeric_time[ND_2_bar == 1], na.rm = TRUE),  # First disaster
    treated2 = ifelse(numeric_time >= first_disaster_year_month, 1, 0)  # Treated after first disaster
  ) %>%
  ungroup()

dataset2_f <- dataset2_f %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    second_disaster_time = nth(numeric_time[ND_2_bar == 1], 2, default = NA)  # Second disaster time
  ) %>%
  ungroup()

# Exclude months after the second disaster
dataset2_filtered <- dataset2_f %>%
  filter(is.na(second_disaster_time) | numeric_time <= second_disaster_time)

test <- output_f %>% filter(year == 2010 | year == 2011) %>% filter(z_value >= 2)
test1 <- test %>% group_by(ADM4_PCODE) %>% summarise(n=n())
test1_b <- test1 %>%  mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", ""))

dataset2_filtered <- dataset2_filtered %>%
  filter(!ADM4_PCODE %in% test1_b$ADM4_PCODE)

main6 <- feols(npa_after1~ treated2 + Mixed_inf2 + Enclave_inf2 +
                 treated2*Mixed_inf2 + treated2*Enclave_inf2
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_filtered)

model_list <- list(main2, main1, main3, main4, main5, main6)
pd <- position_dodge(0.75)

# Extract coefficients and CIs using the broom package
coefficients_df_main <- bind_rows(lapply(model_list, tidy))
coefficients_df_main1 <- coefficients_df_main[c(4,9,14,19,24,29),]

# Create the coefficient plot with CIs using ggplot2
figureA4 <- ggplot(coefficients_df_main1, aes(y = estimate, x = term)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Coefficient estimate (mixed setting)") +
  scale_x_discrete(labels = c("This research", "1.5", "1.7", "1.8", "1.9", "2")) +
  theme_bw(base_size = 14)  +
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6) 

# Extract coefficients and CIs using the broom package
coefficients_df_main2 <- coefficients_df_main[c(5,10,15,20,25,30),]


# Create the coefficient plot with CIs using ggplot2
figureA5 <- ggplot(coefficients_df_main2, aes(y = estimate, x = term)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Coefficient estimate (enclave setting)") +
  scale_x_discrete(labels = c("This research", "1.5", "1.7", "1.8", "1.9", "2")) +
  theme_bw(base_size = 14)  +
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6) 

figure8 <- grid.arrange(
  figureA4, figureA5)
########
########

#### Online Appendix O: A DID model accounting for the treatment phase ####
# First, Convert 'time' to a proper Date object 
dataset1 <- dataset1 %>%
  mutate(year = as.numeric(year),
         numeric_time = (year - 2012))

# Second, code the first hit 
dataset1 <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_year = min(numeric_time[ND_1.6_bar == 1], na.rm = TRUE),  
    # First disaster
    treated = ifelse(numeric_time >= first_disaster_year, 1, 0)) %>%
  # Treated after first disaster 
  ungroup()

# Third, code the second phase  
dataset1_treatment <- dataset1 %>%
  group_by(ADM4_PCODE) %>%
  mutate(
    first_disaster_time = min(numeric_time[ND_1.6_bar == 1], na.rm = TRUE),  # First disaster time
    second_disaster_time = nth(numeric_time[ND_1.6_bar == 1], 2, default = NA),  # Second disaster time
    treatment_phase = case_when(
      numeric_time < first_disaster_time ~ 0,  # Pre-treatment
      numeric_time >= first_disaster_time & (is.na(second_disaster_time) | numeric_time < second_disaster_time) ~ 1,  # Phase 1
      numeric_time >= second_disaster_time ~ 2  # Phase 2
    )) %>%
  ungroup()

# Fourth, I exclude the disaster-affected barangays in 2010 and 2011. 
test <- output_f %>% filter(year == 2010 | year == 2011) %>% 
  filter(z_value >= 1.6) %>% 
  # output_f is the data on the precipitation deviation genereated in 2)
  group_by(ADM4_PCODE) %>% 
  summarise(n=n()) # calculate the affected barangays in 2010 and 2011

dataset1_treatment_filtered <- dataset1_treatment %>%
  filter(!ADM4_PCODE %in% test$ADM4_PCODE)

omod1 <- feols(Tchange_inf2 ~ i(treatment_phase) + Mixed_inf2 + Enclave_inf2 +
                 i(treatment_phase)*Mixed_inf2 + i(treatment_phase)*Enclave_inf2
               | ADM4_PCODE + year, vcov = ~ADM4_PCODE + year,
               data = dataset1_treatment_filtered)

etable(omod1, tex = TRUE)
########

#### Online Appendix Q: Measuring spatial configuration as a continuum ####
# Initialize results list
results <- list()

# First model for `treated*Enemy_around`
model_start <- feols(Tchange_inf2 ~ treated * Enemy_around | ADM4_PCODE + year,
                     vcov = ~ADM4_PCODE + year,
                     data = dataset1_staggeredadoption)

# Extract coefficients and confidence intervals for `treated` and `treated:Enemy_around`
coef_treated <- coef(model_start)["treated"]
coef_interaction <- coef(model_start)["treated:Enemy_around"]
vcov_start <- vcov(model_start)

results[["Enemy_around"]] <- list(
  Coefficient = coef_treated + coef_interaction,
  CI_90_Lower = coef_treated + coef_interaction - 1.645 * sqrt(
    vcov_start["treated", "treated"] +
      vcov_start["treated:Enemy_around", "treated:Enemy_around"] +
      2 * vcov_start["treated", "treated:Enemy_around"]),
  CI_90_Upper = coef_treated + coef_interaction + 1.645 * sqrt(
    vcov_start["treated", "treated"] +
      vcov_start["treated:Enemy_around", "treated:Enemy_around"] +
      2 * vcov_start["treated", "treated:Enemy_around"]),
  CI_95_Lower = coef_treated + coef_interaction - 1.96 * sqrt(
    vcov_start["treated", "treated"] +
      vcov_start["treated:Enemy_around", "treated:Enemy_around"] +
      2 * vcov_start["treated", "treated:Enemy_around"]),
  CI_95_Upper = coef_treated + coef_interaction + 1.96 * sqrt(
    vcov_start["treated", "treated"] +
      vcov_start["treated:Enemy_around", "treated:Enemy_around"] +
      2 * vcov_start["treated", "treated:Enemy_around"]))

# Loop over Enemy_0.x variables
for (i in seq(0.1, 0.9, by = 0.1)) {
  var_name <- paste0("Enemy_", i)
  interaction_term <- paste0("treated:", var_name)
  
  # Dynamically run the model
  formula <- as.formula(paste(
    "Tchange_inf2 ~ treated + Enemy_around +", var_name, "+",
    "treated*Enemy_around + treated*", var_name,
    "| ADM4_PCODE + year"))
  
  model <- feols(formula, vcov = ~ADM4_PCODE + year, data = dataset1_staggeredadoption)
  
  # Extract coefficients
  coef_treated <- coef(model)["treated"]
  coef_enemy_around <- coef(model)["treated:Enemy_around"]
  coef_interaction <- coef(model)[interaction_term]
  
  # Variance-covariance matrix
  vcov_matrix <- vcov(model)
  
  # Calculate variance for the linear combination
  combined_variance <- vcov_matrix["treated", "treated"] +
    vcov_matrix["treated:Enemy_around", "treated:Enemy_around"] +
    vcov_matrix[interaction_term, interaction_term] +
    2 * vcov_matrix["treated", "treated:Enemy_around"] +
    2 * vcov_matrix["treated", interaction_term] +
    2 * vcov_matrix["treated:Enemy_around", interaction_term]
  
  # Store results with both 90% and 95% confidence intervals
  results[[var_name]] <- list(
    Coefficient = coef_treated + coef_enemy_around + coef_interaction,
    CI_90_Lower = coef_treated + coef_enemy_around + coef_interaction - 1.645 * sqrt(combined_variance),
    CI_90_Upper = coef_treated + coef_enemy_around + coef_interaction + 1.645 * sqrt(combined_variance),
    CI_95_Lower = coef_treated + coef_enemy_around + coef_interaction - 1.96 * sqrt(combined_variance),
    CI_95_Upper = coef_treated + coef_enemy_around + coef_interaction + 1.96 * sqrt(combined_variance))}

# Convert results to a data frame
coefficients_df <- data.frame(
  Variable = names(results),
  Coefficient = sapply(results, function(x) x$Coefficient),
  CI_90_Lower = sapply(results, function(x) x$CI_90_Lower),
  CI_90_Upper = sapply(results, function(x) x$CI_90_Upper),
  CI_95_Lower = sapply(results, function(x) x$CI_95_Lower),
  CI_95_Upper = sapply(results, function(x) x$CI_95_Upper))

coefficients_df <- coefficients_df[-11,]
# Adjust the factor levels for the x-axis
coefficients_df <- coefficients_df %>%
  mutate(Variable = factor(Variable, 
                           levels = c("Enemy_around", paste0("Enemy_", seq(0.1, 0.9, by = 0.1))),
                           labels = c(">0", ">10%", ">20%", ">30%", ">40%", ">50%", ">60%", ">70%", ">80%", ">90%")),
         Coefficient = Coefficient * 100,
         CI_90_Lower = CI_90_Lower * 100,
         CI_90_Upper = CI_90_Upper * 100,
         CI_95_Lower = CI_95_Lower * 100,
         CI_95_Upper = CI_95_Upper * 100)

# Plot the coefficients with both confidence intervals
pd <- position_dodge(0.75)

figureA9 <- ggplot(coefficients_df, aes(x = Variable, y = Coefficient)) +
  geom_point() +
  geom_line(group = 1) +
  
  geom_linerange(aes(ymin = CI_95_Lower, ymax = CI_95_Upper),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = CI_90_Lower, ymax = CI_90_Upper),
                 lwd = 1,
                 position = pd) +
  labs(title = "",
    x = "Territories Surrounded by >X% Enemy-Controlled Territories",  # Custom x-axis label
    y = "% Change in a Shift in Territorial Control") +
  theme_minimal()

########

#### Online Appendix R: Co-distribution of two conflict outcomes ####
map <- barangay %>% dplyr::select(ADM4_PCODE, ADM3_PCODE) %>% st_drop_geometry() 
map_1 <- barangay %>% dplyr::select(ADM4_PCODE, ADM3_PCODE) %>% mutate(ADM4_PCODE = str_replace(ADM4_PCODE, "^PH0*", "")) %>% st_drop_geometry() 

datasetr1 <- dataset1 %>% left_join(map)
datasetr2 <- dataset2 %>% mutate(ADM4_PCODE = as.character(ADM4_PCODE)) %>% left_join(map_1)

#aggregate the conflict data into the year level 
datasetr2 <- datasetr2 %>% group_by(ADM3_PCODE, ADM4_PCODE, year) %>%
  summarise(Violent.Conflict_after1 = if_else(sum(npa_after1) > 0, 1, 0)) %>%
  mutate(ADM4_PCODE = paste0("PH", as.character(ADM4_PCODE)))

#merge the territorial shifts data and conflict data
datasetr1 <- datasetr1 %>% filter(ADM3_PCODE %in% unique(datasetr2$ADM3_PCODE)) %>%
  dplyr::select(ADM3_PCODE, ADM4_PCODE, year, Tchange_inf2)

#load the spatial data on administrative unit and merge it with conflict and territorial control shift data 
map <- barangay %>% st_make_valid() %>%
  dplyr::select(ADM3_PCODE, ADM4_PCODE) %>%
  filter(ADM3_PCODE %in% unique(datasetr2$ADM3_PCODE))

dataset_f <- map %>% left_join(datasetr1) %>% left_join(datasetr2)

#Plot the distribution of conflict and territorial shift 
dataset_f_plot <- dataset_f %>% filter(year != 2011) %>%
  group_by(ADM4_PCODE) %>% 
  summarise(Violent.Conflict_after1 = if_else(sum(Violent.Conflict_after1) > 0, 1, 0),
            Control.Change = if_else(sum(Tchange_inf2) > 0, 1, 0)) %>%
  mutate(both = if_else(Control.Change == 1 & Violent.Conflict_after1 == 1, 1, 0))

figure2a <- ggplot() +
  # Plot areas with Control.Change
  geom_sf(data = dataset_f_plot, aes(fill = factor(Control.Change)), size = 0.2, color = "lightgrey") +
  # Add Violent Conflict areas (red)
  geom_sf(data = dataset_f_plot %>% filter(Violent.Conflict_after1 == 1), aes(fill = "Violent Conflict"), color = "red", size = 0.2) +
  # Add areas where both Control.Change and Violent Conflict occur (purple)
  geom_sf(data = dataset_f_plot %>% filter(Control.Change == 1 & Violent.Conflict_after1 == 1), aes(fill = "Both"), color = "purple", size = 0.2) +
  # Define color scales and legend
  scale_fill_manual(
    values = c("grey", "blue", "red", "black"), 
    name = "Legend", 
    labels = c("Neither", "Territorial Control Change", "Battle-related Violence", "Both")
  ) +
  # Apply minimal theme
  theme_minimal() +
  # Set coordinate limits
  coord_sf(xlim = c(125, 126.5), ylim = c(5.2, 8)) +
  # Customize legend position and style
  theme(
    legend.position = c(0.7, 0.25),
    legend.justification = c(0, 1),
    legend.title = element_blank()) +
  # Customize x-axis labels
  scale_x_continuous(breaks = seq(125, 126.5, by = 0.5))
########

#### Online Appendix S: Count models for battle-related violence ####

# Poisson regression with fixed effects and clustered standard errors
dataset2 <- dataset2 %>% arrange(ADM4_PCODE, year, month) %>%
  group_by(ADM4_PCODE) %>%
  mutate(npa_count_1 = dplyr::lead(npa_count, n=1))

smod1 <- fepois(
  npa_count_1 ~ ND_1.6_bar + Mixed_inf2 + Enclave_inf2 +
    ND_1.6_bar*Mixed_inf2 + ND_1.6_bar*Enclave_inf2 |
    ADM4_PCODE + time, vcov = "twoway", data = dataset2)

# Negative Binomial regression with fixed effects
smod2 <- fenegbin(
  npa_count_1 ~ ND_1.6_bar + Mixed_inf2 + Enclave_inf2 +
    ND_1.6_bar*Mixed_inf2 + ND_1.6_bar*Enclave_inf2 |
    ADM4_PCODE + time, vcov = "twoway", data = dataset2)

etable(smod1, smod2, tex = TRUE)
########

#### Online Appendix T: Exploring alternative temporal aggregations of battle-related violence ####
# Figure A10
mlag1 <- feols(npa  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag2 <- feols(npa_after1  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag3 <- feols(npa_after2  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag4 <- feols(npa_after3  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag5 <- feols(npa_after4  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag6 <- feols(npa_after5  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

mlag7 <- feols(npa_after6  ~ treated 
               | ADM4_PCODE + time, vcov = "twoway",
               data = dataset2_staggeredadoption)

model_list <- list(mlag1, mlag2, mlag3, mlag4, mlag5, mlag6, mlag7)

# Extract coefficients and CIs using the broom package
coefficients_df <- bind_rows(lapply(model_list, tidy))
coefficients_df_natural_disasters <- coefficients_df[c(1,2,3,4,5,6,7),]
coefficients_df_natural_disasters$term <-
  c("Violent.Conflict","Violent.Conflict_after1","Violent.Conflict_after2","Violent.Conflict_after3","Violent.Conflict_after4","Violent.Conflict_after5","Violent.Conflict_after6")

# Create the coefficient plot with CIs using ggplot2
pd <- position_dodge(0.75)

figureA10 <- ggplot(coefficients_df_natural_disasters, aes(x = term, y = estimate)) +
  geom_point(position = pd, shape = 19) +
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd)  +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Coefficient estimate of natural disasters",
       color = "") +
  theme_bw(base_size = 14) +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6) +
  coord_flip() +
  scale_x_discrete(labels = c("t", "t+1", "t+2","t+3","t+4","t+5","t+6"))

# Figure A11
#for conditional mode
m4lag1 <- feols(npa ~treated*Mixed_inf2 +
                  treated*Enclave_inf2
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag2 <- feols(npa_after1 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag3 <- feols(npa_after2 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag4 <- feols(npa_after3 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag5 <- feols(npa_after4 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag6 <- feols(npa_after5 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

m4lag7 <- feols(npa_after6 ~ treated*Mixed_inf2 +
                  treated*Enclave_inf2 
                | ADM4_PCODE + time, vcov = "twoway",
                data = dataset2_staggeredadoption)

model_list <- list(m4lag1, m4lag2, m4lag3, m4lag4, m4lag5, m4lag6, m4lag7)

# Extract coefficients and CIs using the broom package
coefficients_df <- bind_rows(lapply(model_list, tidy))

coefficients_df_natural_disasters <- coefficients_df[c(4,9,14,19,24,29,34,5,10,15,20,25,30,35),]
coefficients_df_natural_disasters$term <-
  c("Battle-related Violence","Battle-related Violence_after1","Battle-related Violence_after2","Battle-related Violence_after3","Battle-related Violence_after4","Battle-related Violence_after5","Battle-related Violence_after6","Battle-related Violence","Battle-related Violence_after1","Battle-related Violence_after2","Battle-related Violence_after3","Battle-related Violence_after4","Battle-related Violence_after5","Battle-related Violence_after6")
coefficients_df_natural_disasters$ts <- c("m","m","m","m","m","m","m","e","e","e","e","e","e","e")


# Create the coefficient plot with CIs using ggplot2
figureA11 <- ggplot(coefficients_df_natural_disasters, aes(x = term, y = estimate, color = as.factor(ts))) +
  geom_point(position = pd, shape = 19) +
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error),
                 lwd = 1/2,
                 position = pd)  +
  geom_linerange(aes(ymin = estimate-1.645*std.error, ymax = estimate+1.645*std.error),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  coord_flip() +
  labs(x = "", 
       y = "Coefficient estimate",
       color = "") +
  theme_bw(base_size = 14) +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6,
        legend.position = "bottom") +
  scale_x_discrete(labels = c("t", "t+1", "t+2","t+3","t+4","t+5","t+6")) +
  scale_color_manual(values = c("grey", "black"), labels = c("Interaction of Enclave setting", "Interaction of Mixed setting"))
########

